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Abstract. Previously, we have proposed a direct simulation scheme for colloidal dispersions in a Newto- 
nian solvent [PhysRev.E 71,036707 (2005)]. An improved formulation called the "Smoothed Profile (SP) 
method" is presented here in which simultaneous time-marching is used for the host fluid and colloids. The 
SP method is a direct numerical simulation of particulate flows and provides a coupling scheme between the 
continuum fluid dynamics and rigid-body dynamics through utilization of a smoothed profile for the col- 
loidal particles. Moreover, the improved formulation includes an extension to incorporate multi-component 
fluids, allowing systems such as charged colloids in electrolyte solutions to be studied. The dynamics of the 
colloidal dispersions are solved with the same computational cost as required for solving non-particulate 
flows. Numerical results which assess the hydrodynamic interactions of colloidal dispersions are presented 
to validate the SP method. The SP method is not restricted to particular constitutive models of the host 
fluids and can hence be applied to colloidal dispersions in complex fluids. 

PACS. 47.11.-j Computational methods in fluid dynamics - 82.70.-y Disperse systems; complex fluids - 
82.20.Wt Computational modeling; simulation 

1 Introduction 

Interparticle interactions in colloidal dispersions mainly 
a e-mail: ynakayama@chem-eng.kyushu-u.ac.jp consist of thermodynamic potential interactions and hy- 
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drodynamic interactions [UIllE] . Whereas the former in- ticle and the boundary value problem of the solvent con- 

teractions occur in both static and dynamic situations, tinuum equations was solved based on the irregular mesh, 

the latter occur solely in dynamic situations. Although For applying this scheme to many-body system, inaccessi- 

thermodynamic interactions have been studied extensively ble computational resources is inevitable. In Refs. [T6jJ17t 

and summarized as a concept of the effective interaction [4], H51TT9] . they adopted the lattice Boltzmann (LB) method 

the nature of dynamic interactions is poorly understood, for explicitly solving the solvent hydrodynamics. Although 

Since the hydrodynamic interaction is essentially a long- the LB method has many advantages for solving large sys- 

range, many-body effect, it is extremely difficult to study terns, compared to the standard discretization schemes, 

its role using analytical methods alone. Numerical simu- the applicability to various constitutive equations for the 

lations can aid the investigation of the fundamental role complex fluids other from the Newtonian fluids is unveiled, 

of hydrodynamic interactions in colloidal dynamics. Concerning the interaction between colloidal particles and 

the solvent, authors in Refs [TBlll7j . arranged finite inter- 



In recent decades, various simulations for particulate 
flow have been developed for most simple situations in 
which the host fluid is Newtonian [5,6,7,8,9,10, lT|J 12^1X3] . 
However, these schemes can not necessarily be applied to 
problems with non-Newtonian host fluids or solvents with 
internal microstructures, which are practically more im- 
portant cases. Hydrodynamic simulations of ions and/or 
charged colloids have been proposed by making use of In this article, we propose a direct simulation scheme 



action points on the surface of the colloids on which a fric- 
tional coupling was assumed. In this coupling scheme, fric- 
tion parameters and the number of the interaction points 
were determined phcnomcnologically and theoretical basis 
to determine the parameters in general systems was not 
clearly stated. 



some of the above schemes 14,15,16,17,18,19 . Nonethe- for colloidal dispersions which is applicable to most consti- 

less, the tractability and/or physical validity of their mod- tutive models of the host fluids. We call it the "Smoothed 

eling remain controversial. In Ref. [14j . in solving the mo- Profile (SP) method" since the original sharp interface 

tion of the ionic solutes, the solvent hydrodynamic interac- between the colloids and solvent is replaced with an effec- 

tion was incorporated by Rotne-Pragar-Yamakawa-type tive smoothed interface with a finite thickness [71H2U]. We 

mobility tensor, which accounted for a long-range part of formulate a computational method to couple the particle 

pair interaction. This scheme misses many-body and/or dynamics and hydrodynamics of the solvent. A fixed grid 

near-field hydrodynamic interactions between the solutes is used in both the solvent and the particle domains. Intro- 

and surfaces, which can have an effects on the behavior duction of a smoothed profile makes it possible to realize 

of ions near surface. In Ref. [15], the computational mesh stable and efficient implementation of our scheme. The nu- 

was arranged to express round shape of a colloidal par- merical implementation for Newtonian solvents and elec- 
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trolyte solutions as specific examples of multi-component Here, for the sake of simplicity, we do not deal with the 
fluids is outlined. Various test cases which verify the SP cross diffusion of different solutes. The conservation of mo- 
method and assess the hydrodynamic interactions are pre- mentum implies the velocity of solvent follows the Navier- 
sented. Stokes equation of incompressible flow with the source 



2 Dynamics of a multi-component solvent 
and colloids 



term from solutes: 

V • v = 0, (4) 



p (d t + v ■ V) v = -Vp + tiV 2 v - C a Vfi a + V • s, 

a 

2.1 Hydrodynamics of multi-component fluids (5) 



where p is the total mass density of the fluid, p is the pres- 
sure, r\ is the shear viscosity of the fluid, and s is a random 
stress satisfying the fluctuation-dissipation relation [2Tj : 



We first give a brief description of multi-component fluid 
equations, looking at electrolyte solutions as a specific ap- 
plication. Consider N (possibly ionic) solute species that 
satisfy the law of conservation for every concentration, C a , {sn,(x, t)sji(x', t ))& = 2kBTrj(8ij8ki + 5u6kj)8(x — x)8(t — t). 
of the ath species: (6) 

The above set of equations is closed when a set of chem- 

d t C a + V • C a v a + V • g a = 0, (1) 

ical potentials {p a } is given, and describes the dynamics 

where v n is the velocity of the ath solute and q„ is a ran- „ ,,. ,n-in,i -n i. r 

a of a multi-component fluid, ior the specific application of 

dom current. Since the inertial time scales of the solute , , , , , , . . , ,, „ . , T 

an electrolyte solution, we consider the Foisson-JNernst- 

moleculcs are extremely small, the velocity of the ath so- „, . , . „ , , , . , , . , 

Planck equation for the chemical potential: 

lute can be decomposed into the velocity of the solvent v 
and the diffusive current arising from the chemical poten- 
tial gradient V ' p a as: 



p a {{C u . . . , C N }) = k B TlogC a + Z a e($ - E ext ■ x), 

(7) 

eV 2 $ = -p e , (8) 

This equation describes the Poisson-Boltzmann distribu- 
tion for ions at equilibrium, where Z a is the valence of the 



v a =v- r a Vp, a , (2) 

where r a ksT is the diffusivity of the ath ion, Icb is the 
Boltzmann constant and T is the temperature. The ran- 

ath ion, e is the elementary charge, <P is the electrostatic 
dom current should satisfy the following fluctuation-dissipation 

potential, E ext is the external field, e is the dielectric con- 
relation [2"T] : 

stant of the fluid, and p e is the charge density field. This 
(9a,i(x,t)g^j(x' ,t')) = 2(/csT) 2 r a 8 a p8ij8{x — x')S(t — t'). get of equations corresponds to the electrokinetic equa- 
ls) tions which appear in standard textbooks [3]. 
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2.2 Colloids in electrolyte solutions eluded in the term F\ x , and the external field on the 



The colloid dynamics are maintained by the force exerted 
by the solvent. Consider monodisperse spherical colloids 
with a radius a, a mass M p , and a moment of inertia 
I p . Momentum conservation between the fluid and the 
zth colloid implies the following hydrodynamic force and 
torque: 

Ff = J (dS< ■ a), Nf = J (x - Rt) x (dS< ■ <x), (9) 



colloids is accounted for by Ff . 



Having colloidal particles with a finite volume pro- 
vides the relevant boundary conditions to the hydrody- 
namic equations. For the solvent velocity, a no-slip con- 
dition is assigned, such that v = V i + fli x n with 
Ti = x Ri for the zth colloid. For the concentration field, 
a no-penetration condition is assumed, giving n • V^ Q = 0, 
where n represent the unit normal to the surface of the 



where R t is the center of mass, / dS t (. . . ) indicates the colloids. Coupling the hydrodynamics of the solvent with 
surface integral over the ith colloid, and a is the stress thc d y nam ics of the colloids defines the moving-boundary- 
tensor of thc fluid. In terms of thc clcctrokinctic equations, condition problem above. The usual numerical techniques 
the stress reads as °^ P ar ^ial differential equations arc hopeless in dealing 

with the dynamical evolution of many colloids since the 

a = -pi + a' + cr st + s, (10) 

sharp interface at the surface of the colloids moves and 
in which a' = r\ (\7v + (Vt>) T ) is thc dissipative stress, henceforth the mesh points at which the boundary con- 
and er st = e ^EE - (\E\ 2 /2)7 j is the Maxwell stress for dition is assigned vary with each discrete time step. Thc 
thc electric field E = -\/<P+E ext , where I is thc unit ten- moving boundary condition leads to huge computational 
sor. The evolution of the colloids follows Newton's equa- costs. 

* lons ' In contrast, the SP method formulates an efficient scheme 

• . . for this kind of moving-boundary-condition problem, and 

Ri = V i, (11) 

incurs the same level of computational cost as is required 

M p V t = Ff + F c l +F e l xt , (12) 

for solving a uniform fluid. The typical computational 

I p -f2i = N? + Nf xt , (13) 

costs for the fluid and particles are assumed to scale to 
where F ext and N ext are the external force and torque, their degrees of freedom. In thc SP method, regular mesh, 
respectively, and Fl is the force arising from the core po- but not body-fitted mesh, can be used, which indicates 
tential of the particles, which prevents the colloids over- that the inclusion of the dispersed particle phase does not 
lapping. Hereinafter, the soft-core potential of the trun- induce the increase of the grids. In (i-dimensional system, 
cated Lennard- Jones potential is adopted for F\. For the we assume that discretized space contains N d grids and 
charged-colloid system specifically, the buoyancy is in- the variables of the fluid phase are linked to thc grids. 
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N p particles are dispersed in the system and their volume tration field, is forced to be zero through multiplication 
scale to N p a d . The number of the particles are at most by (1 — 0). 

N p a d < (NAx) d where Ax is the lattice spacing, which re- The advection of 4> is solved via Ri — Vj and by map- 
suits in N p /N d < (Ax/a) d . This inequality indicates that ping . . . , Rn p } to </>. Henceforth, the volume of the 
the computational costs for the dispersed phase tracking fluid and/or the solid is strictly conserved and no numcr- 
is at largest (Ax/a) d times for the fluids. In typical case ical diffusion of occurs. In the SP method, the funda- 
of a I Ax = 5, the fraction of the computational costs for mental field variables to be solved are taken as the total 
the particles roughly estimated less than several percent, velocity v rather than Vf, and C* rather than C a . This 
thus most of the computation should be for solving the choice of variables yields great benefits in terms of allowing 
uniform fluid. efficient and stable time evolution. The evolution equation 

of v is derived based on momentum conservation between 
3 Computational algorithm *he ^ m< ^ anc ^ * ne particles, and the rigidity of the parti- 

cle velocity field v p . For the solute concentration, C a and 

In the SP method, quantities are defined over the entire 

C* differ in terms of whether they exhibit abrupt varia- 

domain, which consists of the fluid domain and the particle 

tion at the interface of the colloids or not: especially for 

domain. To designate the particle domain, we introduce 

£ — > 0, C a has discontinuity at the solid-fluid interface, 

a concentration field for the colloids, given as (f>(x,t) = 

but C* should not. Since C a exhibits abrupt variation in 

Hi^i <t>i{ x it)i where fa E [0, 1] is the profile field of the 

the spacial scale of £, in order to solve its advection, it is 

ith particle, which is unity at the particle domain, zero 

necessary to stabilize the evolution over time and to pre- 

at the fluid domain, and which has a continuous diffuse 

vent numerical diffusion and penetration of C a into the 

interface of finite thickness £ at the interface domain. With 

particle domain. Compared with C a , auxiliary concentra- 

the field 0, the total velocity field and concentration fields 

tion C* is independent of the abrupt variation arising from 

of the solutes are denned as 

(1 — <f>), and finite value of C* in the particle domain is 
v = (1 — 4>)vf + (j)v p , (14) allowed. Therefore, it is numerically much easier to solve 

C a = (1 — <f>)C*, (15) the advection equation of C* than C a . 

where (1 — <p)v f represents the velocity field of the fluid, 

3.1 Discretization in time 

and <j>Vp(x,t) = E^i [Vi(t) + fli(i) x n(t)] is 

the velocity field of the colloids. The auxiliary concentra- The time-discretized evolution equations are derived as 
tion field C* is introduced, which can have a finite value follows. To simplify the presentation, we neglect random 
in the particle domain, whereas C a , the physical concen- currents. As an initial condition at the nth discretized time 
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step, the position, velocity, and angular velocity of the with the charge density field, = J2 a ^aeC*' n+ 



colloids, {!?", VI, fl 1 }} (i = 1, . . . , N p ), are mapped to cj) n 
and <p n Vp and the following conditions are set: v n = (1 — 
(j) n )vf + 4> n Vp, satisfying the incompressibility condition 
on the entire domain, V • v n = 0, and C*' n , satisfying the 
charge neutrality condition, J dx(l — <p n ) J2 a Za^C" 1 + 
J dx |V0| ea e = 0, where V</>i| ea e represents the surface 
charge distribution of the colloids. The current for the 
auxiliary concentration field is defined as 



C* a v a = C* a v + nn) ■ C* {-r a V^ a ) . 



where n(x,t) is the unit surface-normal vector field which 
is defined on the interface domain with a finite thickness £. 
In this definition of the current, the no-penetration condi- 
tion is directly assigned. The auxiliary concentrations are 
advected by this current as 



a 



*,n+l 



c: 



rt n +h 

J ds\7-C* a v a , (17) 



where h is a time increment, and t n = nh is the nth 
discretized time. The total velocity field is updated using 
a fractional step approach. First, the advection and the 
hydrodynamic viscous stress are solved, 

1 



rt n +h 

V* =v n + J dsV 
rt n +h 



(—pi + <j') — vv 



dsVt, 



(18) 
(19) 



with the incompressibility condition, V ■ v* =0. Along 
with the advection of the total velocity, the particle po- 
sition is updated using the particle velocity. The electro- 
static potential for the updated particle configuration is 
determined by solving the following Poisson equation, 

V 2$n+1 = _ p «+l/ £; (20) 



| V0™ +1 1 ea e . The momentum change as a result of the 
electrostatic field is solved as 



v** = v* - hp r ^ +1 V<P n+1 



(21) 



At this point, the momentum conservation is entirely solved 
for the total velocity field. The rest of the updating pro- 
cedure applies to the rigidity constraint on the particle 
(16) velocity field. 



The hydrodynamic force and torque on the colloids 
exerted by the fluid are derived by assuming momentum 
conservation between the colloids and the fluid. The time- 
integrated hydrodynamic force and torque over a period 
h are equal to the momentum change over the particle 
domain: 



rtn+h 

^ dsFf(s) 

rt n +h 

jf dsN?( S ) 



J dxp^ +1 (v** - «») , (22) 

J dx [r- +1 x p^ +1 („** - v;)] . 

(23) 



With this and other forces on the colloids, the particle 
velocity and angular velocity are updated as 



yn+l _yn, M ~l 



n n+1 = n n + 1 1 ■ 



tn+h 



dsF 



H 



+M- i jf* B i(Jt+Ff rt ) 

(24) 



rtn+h 




I dsiYf 







rtn+h 



(25) 
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The resultant particle velocity field n+1 u™ +1 is directly is used for the velocity over the particle domain. For corn- 
enforced on the total velocity field as parison, in the Fluid Particle Dynamics (FPD) method [7], 



v n+1 = v* 



ds(f>f p 



the large viscosity is introduced for the fluid particle rj c 



t„+h 

dscf>f p 



h n+l / n+1 



(3> 77) to realize the rigidity constraint. This means that 
v **) - ~ V Pp' ( 27 ) the required time increment for FPD should be very small, 



P 

where <j>f p represents the force density field which im- Le -> f/»fc(« 1), as compared with that used in the SP 
pose the rigidity constraint on the total velocity field, method. 
The pressure due to the rigidity of the particle is deter- 
mined by the incompressibility condition, V • v n+1 = 0, 
which leads to the following Poisson equation for p p , viz, 
V 2 p p = £ V • [<f) n+1 (v n p +1 - v**)] . We note again that, on 
the l.h.s. of Eqs.fl^,©, and ([27]), the integrands Ff 



N? , and (j)f„ are not explicitly calculated but their time . . , 

F condition m the advection-diftusion equation of the so 



A similar discussion on the restriction of the no-penetration 

con 

integrals are solved. In other words, the solid-fluid inter- . 

lutes can be outlined. One of the simplest treatments of 

actions are treated in the form of the momentum change, . ..... ...... 

the no-penetration condition m the particle domain is the 

namely, momentum impulse 



penalty method adopted in refs [23 , 24J , in which an arti- 
ficially large potential barrier is introduced for the solutes 

3.2 Restriction on a time increment 

in the particle domain. The artificial potential should at 
To enable spatial discretization of the hydrodynamic equa- least be larger than the other chemical potentials in order 
tions, any standard scheme, such as the finite difference to realize no-penetration of the solutes. This restriction 
method, finite volume method, finite element method, spec- means that the artificial potential requires a smaller h. 
tral method, lattice Boltzmann discretization and so forth, Although this strategy is physically consistent, it is nu- 
can be used. The SP method basically defines a coupling merically inefficient. In contrast, in the SP method the 
scheme between the hydrodynamic equations for the sol- advection-diffusion requires no additional time scale for 
vent and the equations for the discrete colloids. Since the the inclusion of colloids since the no-penetration condition 
treatment of the rigidity constraint of the particle velocity for the solutes is directly assigned to the solute current in 
does not introduce an additional time scale, the restriction the finite interface domain. From the above discussion, we 
to a time increment h is the same as that in uniform fluid can see that the SP method provides us with much higher 
cases. This is advantageous as compared to the methods numerical efficiency than other methods proposed for di- 
adopted in refs [71l22j. where a large viscosity or elasticity rect numerical simulation of colloidal dispersions. 
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Fig. 1. Snapshot of a sedimenting colloid of radius a = kAx. 
The arrows indicate the velocity field. 



simulation 

Numerical (Zick & Homsy 1 982) 
low-<p asymptotics (Hasimoto 1 959) 




0.2 0.3 0.4 0.5 
volume fraction 



Fig. 2. Drag coefficient of a periodic array of spheres in steady 
Stokes flow as a function of volume fraction ip solved using the 
SP method as compared with the analytic result |25j (Solid 
line) and low-y asymptotics (dashed line). In this case, 
£ is set to unity in the lattice unit. The range of the volume 
fraction ip was obtained by changing the radius a from 4 to 
15.5 for a lattice size Ax with L = 32Ax fixed. 

4 Results and Discussions 

4.1 Stokes drag on a periodic array of spheres in a 
Newtonian fluid 

To validate the SP method quantitatively, we measured 
the steady state drag force on a periodic array of spheres 



in a Newtonian solvent. The velocity distribution around 
the colloid is depicted in FigfT] In general, flow around 
a colloid occurs as creep-flow with a Reynolds number of 
Re = aV/v <C 1. Figure [2] shows the drag coefficient Q((p), 
defined as 



F. 



D 



QitrjaV 



(28) 



Q(<p) ' 

where if = (4/3)7r(a/L) 3 is the volume fraction in a cubic 
box of volume L 3 . The effect of the boundary condition 
extends very far in creep-flow, as the higher the volume 
fraction, the larger the drag. Comparison of the drag co- 
efficient given by the SP method with the analytic result 
from the Stokes equation by Zick and Homsy |25j verifies 
the validity of the SP method over the entire range of tp. 

4.2 Lubrication interaction in a finite system 

One of the most important effects of solvent flow is the 
lubrication interaction between nearby particles with rela- 
tive motions. The exact solution of the Stokes equation for 
two isolated spheres has been found [27] and has provided 
much insight into the basic physics of colloidal suspen- 
sions. However, its application to many-particle systems 
through the method of pairwise addition requires care. 
For quantitative prediction of the rheology of concentrated 
suspensions, numerical results have identified many differ- 
ences, whether the shear mode of the lubrication interac- 
tion is included or not 28 , 29J . There exists a fundamental 
ambiguity in the application of the analytic expression of 
two isolated spheres to a many-particle system in the finite 
domain using pairwise addition. 
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Fig. 3. Snapshot of two approaching colloids of radius a = 
4 Ax. The arrows indicate the velocity field. The color Map 
around the colloids represents the pressure distribution: a 
change in color from red to blue corresponds to a change from 
high to low pressure. 



We computed the squeeze lubrication interaction be- 
tween two approaching spheres in a finite system. The ve- 
locity and pressure distributions are depicted in Fig[3] Fig- 
ure ID^a) shows the normalized approach velocity of a pair 
of particles versus the gap h between two equal spheres 
solved by the SP method as compared with other theo- 
ries. The two asymptotic solutions at h -C a and h^$> a are 
from the exact solution of the isolated pair [57] and Rotne- 
Prager-Yamakawa tensor [30] , respectively. The Stokesian 
Dynamics (SD) solution [5] is based on an interpolation 
of these two asymptotic solutions. The simulation result 
nicely reproduces not only the two asymptotic regimes but 
also the crossover which occurs at h/a ~ 0.7. It is found 
that SD underestimates the mobility in this intermediate 
regime. This underestimation is due to the approximation 
adapted in SD in which the components are just asymp- 
totic two-body solutions. The result confirms the relevance 




simulation 
RPY tensor(Ewald summed) - 
2 body exact resistance-tensor R 2B 
Stokesian Dynamics ■ 



0.0 0.5 1.0 1.5 2.0 2.5 3.0 

h/a 



(hi T 1 
(b) y 2 

:isolated single sphere ■ 



0.1 0.2 0.3 0.4 

volume fraction 



Fig. 4. (a) Relative velocity between two approaching spheres 
versus the gap between the sphere surfaces (symbols). The 
slight oscillations in the SP results are a result of the finite 
lattice spacing. In this case, £ is set to unity in the lattice 
unit. Theoretical curves are shown for far-field asymptotics us- 
ing the Rotne-Prager-Yamakawa tensor (dotted line) [30] . a 
near-field exact solution (dashed line) [27] . and Stokesian Dy- 
namics (solid line) [5] . (b) The two coefficients representing the 
squeeze interaction (scaled to those at infinite dilution) versus 
volume fraction. 71 (square) represents the one-body drag and 
72 (circle) represents the two-body squeeze interaction due to 
relative motion. 

of our simulation, thus demonstrating the importance of 
the hydrodynamic interaction in a finite system. 

We discuss the dependence of the lubrication inter- 
action on the volume fraction. By assuming the func- 
tional form \Vi-V 2 \/F = 1/ (A (if) + 2P 2 O) /h) from 
lubrication theory, where = 6^77071 (ip) is the one- 



10 Yasuya Nakayama et al.: Simulating (electro) hydrodynamic effects in colloidal dispersions: smoothed profile method 

body drag coefficient and ^((p) = (3nria 2 /2)j2( l fi) is the The skewed ion-distribution gives rise to polarization in 
squeeze coefficient, the effect of the volume fraction is rep- the electric double-layer. Moreover, double-layer deforma- 
resented by 71 (y) and 72 (<p). These friction coefficients tion induces an electro-osmotic flow, which makes the flow 
can be extracted by fitting the curve of Fig|4ja) for each ip. around a charged colloid different from that of a neutral 
The reduced coefficients 71 and 72 are plotted in Figj4f b) colloid, 
as a function of the volume fraction. The solid line in 
FigHIb) is from FigO Although a periodic array of spheres 

In this simulation, a periodic array of charged col- 
exhibits different flow geometry from the case of two ap- 

loids in a 1:1 electrolyte solution under gravity was com- 

proaching spheres, the volume fraction dependence of 71 

puted. We specify the valence of the colloids first, and the 

of the two cases is comparable. Moreover, the squeeze co- 

counterion in the host fluid assures the charge neutral- 
efficient 72 is found to be a decreasing function of the vol- 

ity of the whole system. The bulk concentrations of the 

ume fraction. In other words, the squeeze mode is most 

1:1 electrolyte C± were determined by specifying the De- 
enhanced at infinite dilution. In the literature [2], it is . 

bye length k" 1 = {Att\ b {C + + C_)}~ 1 . The Bjerrum 

pointed out that the squeeze coefficient is at least smaller 

length 4ttXb = e 2 /efcgT was set to unity and £ was cho- 

than that of the exact solution for isolated pairs. These 

sen as twice the lattice unit to resolve the surface charge 

results further validate the SP method. 

distribution, which is represented by |V0|. For simplicity, 

Although the SP method itself is an efficient scheme as 

the counterion and coion were set to the Schmidt number 

a direct simulation of utility in constructing a more coarse- 

0.5, which choice did not affect the qualitative aspect of 

grained model for suspensions, such as in dissipative parti- 

the following results but only control the transient time 

cle dynamics, constitutive modeling, etc., the calculation 

to the steady state. The Debye length was chosen so that 

by direct simulation also gives fundamental information 

the effective radius of the double-layer a + n^ 1 was smaller 

about the hydrodynamic interactions. 

than half of the system size L/2. In order to clarify the 
situation in terms of a typical colloid science parameter, 

4.3 Sedimenting charged colloids 

nondimensional zeta potential eQ/ksT is shown in Fig0 
As an example of the specific application of the method to which should be determined by specifying na and the va- 
a multi-component fluid, we compute the hydrodynamic lence of the colloid Z . The corresponding dimension of 
drag for sedimenting charged colloids in the absence of the zeta potential C for this simulation at 20° C is of the 
an external electric field. In this case, the sedimenting order of lOmV. We computed the linear response regime 
charged colloids induce a flow that determines a charge for gravitational driving to observe the effect of electro- 
distribution which differs from the equilibrium case [51131] . osmosis on hydrodynamic drag. 
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We plot the sedimentation velocity scaled by that of a cays as v oc e~ Kr /r [35] -These factors account for why the 

neutral colloid as a function of the inverse Debye length flow patterns in Figs. ([I]) and ([6]) resemble one another. In 

scaled by the radius of colloid in Fig. [5] As the valence other words, the electrohydrodynamic interaction should 

of the colloid or the zeta potential increases, the sedimen- be a pronounced effect in small finite systems as is seen in 

tation velocity reduces to that of the neutral colloid. The neutral systems. 

hydrodynamic drag of the electrolyte solution was most We note that, at a finite volume fraction and in the 

enhanced when the Debye length was comparable to the ran ge of the Debye length, the charged colloid dynamics 

size of the colloids. These facts qualitatively agree with cou ld be effectively assessed through the SP method. Fur- 

the analytic result at infinite dilution [32[33] . thcr application of the SP method to electrophoresis of 

The charge-density and the velocity distributions when concentrated suspensions is reported elsewhere [3S] . 
Ka = 1 and Z = 1000 are depicted in Fig(7j In this case, 

3.5 

the charge distribution was largely uniform and thus the 30 

2 5 

counteracting electrostatic effect of the double-layer polar- 
ization does not have much effect. Therefore, the enhanced § ( g 

electrohydrodynamic drag can mainly be attributed to the 

0.5 

friction between the solvent and the ions. Because of this 
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electrohydrodynamic coupling of the transfer in momen- 
tum, the solvent-flow pattern around the colloid was mod- Fi S- 5 - Non-dimensionalized zeta potential as a function of 

ifred from the case of a neutral colloid. As a result, the vis- the inverse Deb y e len S th Ka and the valence of the colloid Z 

j n -j i j mi • l computed using the semi-analytic formula of Ohshima-Healy- 

cous drag on the colloid was enhanced. 1ms mechanism ° J 

r , , . iii . n n White 1371 . The dotted lines are guides to the eye. 

of enhanced viscous drag by electro-osmotic flow generally 

exists in colloids in electrolyte solutions. 

It has been known that in an infinitely dilute system 

(or the thin double-layer limit) the velocity decays as r~ 3 5 Conclusions 

in the region of sr > 1 in contrast to the r -1 decay of 

an infinitely dilute neutral system, where r is the radial We have presented a new simulation scheme for colloidal 

distance from the center of the colloid |34U35j . However, dispersions in a solvent of a multi-component fluid, which 

in the system size adapted in our simulation this asymp- we call the "Smoothed Profile (SP) method" . The SP 

totic regime was not reached. For nr < 1, we have the method improves and extends a previous method proposed 

screened hydrodynamics regime, where the velocity de- by the authors 38,20,39 . 
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Fig. 6. Sedimentation velocity of a periodic array of colloids 
of different valences Z in electrolyte solutions as a function 
of the inverse Debye length na. The ordinate is reduced by 
the sedimentation velocity of neutral colloids at the same vol- 
ume fraction of 0.008 (a = AAx, L = 32Ar). The effective 
volume fraction including the electric double-layer defined by 
(4/3)7r{(a + «: _1 )/L} 3 was chosen to be less than unity. Dotted 
lines are guides to the eye. 




The description of the colloidal systems is based on 
the Navier-Stokes equation for the momentum evolution 
of the solvent flow, the advection-diffusion equation for 
the solute distribution, the rigid body description of col- 
loidal particles, with dynamical coupling of all these el- 
ements. Based on momentum conservation between the 
continuum solvent fluid and the discrete rigid colloids, 
the time-integrated hydrodynamic force and torque are 
derived. This expression of the mechanical coupling be- 
tween a fluid and particles is well-suited for numerical sim- 
ulations in which differential equations are discretized in 
time. With this formulation relating the solid-fluid inter- 
action, standard discretization schemes for uniform fluids 
are utilized as is; no special care is needed for the solid- 
fluid boundary mesh. Since the hydrodynamic interaction 
is solved through direct simulation of the solvent fluid, the 
many-body effect can be fully resolved. 



The utility of the SP method was assessed in various 
test problems, not only using simple fluids but also simu- 
lating charged colloids in electrolyte solutions. The results 
confirmed that the SP method is effective for studying 
Fig. 7. Snapshot of a sedimenting charged colloid when Ka = 1 the dynamical behavior of colloids. Although we have fo- 
and Z = 1000. The arrows indicate the velocity field. The color cuged Qn systems of spher ical colloids in simple fluids and 

map around the particle represents the charge density field: a „ . AT , n1 iij.ij.c-ti- n u 
r Poisson-JNernst-Planck electrolytes (i.e. Poisson-Boltzman 



change in color from red to blue corresponds to a change from 
high to low charge density. 



level description of electrolytes) , application to other types 
of macromolecules of other shapes, such as disks [T5], rods, 
and others, and other constitutive models for solvents is 
straightforward. 



Yasuya Nakayama et al.: Simulating (electro) hydrodynamic effects in colloidal dispersions: smoothed profile method 13 
References 16. V. Lobaskin, B. Diinweg, C. Holm, J. Phys.: Condens. 



Matter 16, S4063 (2004) 

17. A. Chatterji, J. Horbach, J. Chem. Phys. 122, 184903 
(2005) 

18. F. Capuani, I. Pagonabarraga, D. Frenkel, J. Chem. Phys. 



1. J. Happel, H. Brenner, Low Reynolds number hydrodynam- 
ics: with special applications to particulate media, 2nd edn 
(Martinus Nijhoff, Dordrecht, 1983) 

2. S. Kim, S.J. Karrila, Microhydrodynamics : principle, 
and selected applications (Butterworth-Heinemann, Lon- 121, 973 (2004) 

don 1991) ^' Capuani, I. Pagonabarraga, D. Frenkel, J. Chem. Phys. 

3. W.B. Russel, D.A. Saville, W.R. Schowalter, Colloidal Dis- 124 > 124903 (2006) 

Persians (Cambridge University Press, Cambridge, Eng- 20 - Y - Nakayama, R. Yamamoto, Phys. Rev. E 71, 036707 

land, 1989) ( 2005 ) 

4. C.N. Likos, Phys. Rep. 348, 267 (2001) 21. L.D. Landau, E.M. Lifshitz, Fluid mechanics (Pergamon 

5. J.F. Brady, G. Bossis, Annu. Rev. Fluid Mech. 20, 111 Press, London, 1959) 

(1988) 22. C.S. Peskin, D.M. McQueen, J. Comput. Phys. 81, 372 

6. A. Malevanets, R. Kapral, J. Chem. Phys. 110, 8605 (1989) 

( 1999 ) 23. J. Dzubiella, H. Lowen, C.N. Likos, Phys. Rev. Lett. 91, 

7. H. Tanaka, T. Araki, Phys. Rev. Lett. 85, 1338 (2000) 248301 (2003) 

8. T. Kajishima, S. Takiguchi, H. Hamasaki, Y. Miyake, ^ R R Takeshit£l) T Araki] H Tfmaka; j phyg 
JSME Int. J., Ser. B 44(4), 526 (2001) Condens Matto ^ ^ {2QM) 

9. H.H. Hu, NA. Patankar, MY. Zhu, J. Comput. Phys. 192, 

25. AA. Zick, CM. Homsy, J. Fluid Mech. 115, 13 (1982) 

26. H. Hasimoto, J. Fluid Mech. 5, 317 (1959) 

27. D.J. Jeffrey, Y. Onishi, J. Fluid Mech. 139, 261 (1984) 

28. R.C. Ball, J.R. Melrose, Physica A 247, 444 (1997) 

29. J.R. Melrose, R.C. Ball, J. Rheol. 48, 937 (2004) 

30. C.W. Beenakker, J. Chem. Phys. 85, 1581 (1986) 

31. R.F. Probstein, Physicochemical Hydrodynamics : An In- 
troduction, 2nd edn. (John Wiley & Sons, New York, 2003) 



427 (2001) 

10. R. Glowinski, T.W. Pan, T.I. Hesla, D.D. Joseph, 
J. Periaux, J. Comput. Phys. 192, 363 (2001) 

11. A.J.C. Ladd, R. Verberg, J. Stat. Phys. 104, 1191 (2001) 

12. J.T. Padding, AA. Louis, Phys. Rev. Lett. 93, 220601 
(2004) 

13. M.E. Cates, K. Stratford, R. Adhikari, P. Stansell, J.C. De- 
splat, I. Pagonabarraga, A.J. Wagner, J. Phys.: Condens. 
Matter 16, S3903 (2004) 32 ' F " Booth ' J ' Chem ' ^ 22 ' 1956 ( 1954 ) 

14. Y.W. Kim, R.R. Netz, Europhys. Lett. 72, 837 (2005) 33 - H - Ohshima, T.W. Healy, L.R. White, R.W. O'Brien, J 

15. T. Yamaue, M. Sasaki, T. Taniguchi, Multi-Phase Dy- Chem. Soc, Faraday Trans. 2 80, 1299 (1984) 
namics Program "Muffin" User's Manual, |http://octa.jp| 34. J.L. Anderson, Annu. Rev. Fluid Mech. 21, 61 (1989) 
(2005) 35. D. Long, A. Ajdari, Eur. Phys. J. E 4, 29 (2001) 



14 Yasuya Nakayama et al.: Simulating (electro) hydrodynamic effects in colloidal dispersions: smoothed profile method 

36. K. Kim, Y. Nakayama, R. Yamamoto, Phys. Rev. Lett. 
96, 208302 (2006) 

37. H. Ohshima, T.W. Healy, L.R. White, J. Colloid Interface 
Sci. 90, 17 (1982) 

38. R. Yamamoto, Phys. Rev. Lett. 87, 075502 (2001) 

39. K. Kim, R. Yamamoto, Macromol. Theory Simul. 14, 278 
(2005) 



